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Abstract 

We present a study of the multiplicities, of the lateral distributions 
and of the ratio of the electromagnetic to the hadronic components in the 
air showers, generated by the collision in the atmosphere of an incoming 
high energy cosmic ray and mediated by the formation of a mini black 
hole, predicted in TeV scale gravity models with large extra dimensions. 
The analysis is performed via a large scale simulation of the resulting cas- 
cades over the entire range (10^^ - lO^^eV) of ultra high initial energies, 
for several values of the number of large extra dimensions, for a variety of 
altitudes of the initial interaction and with the energy losses in the bulk 
taken into account. The results are compared with a representative of 
the standard events, namely the shower due to the collision of a primary 
proton with a nucleon in the atmosphere. Both the multiplicities and the 
lateral distribution of the showers show important differences between the 
two cases and, consequently, may be useful for the observational charac- 
terization of the events. The electromagnetic/hadronic ratio is strongly 
fluctuating and, thus, less decisive for the altitudes considered. 
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1 Introduction 



In the almost structureless fast falling with energy inclusive cosmic ray spectrum, 
two kinematic regions have drawn considerable attention for a long time [1]. 
These regions are the only ones in which the spectral index of the cosmic ray 
flux shows a sharper variation as a function of energy, probably signaling some 
"new physics" , according to many. These two regions, termed the knee and the 
ankle [2] have been puzzling theorists and experimentalists alike and no clear and 
widely accepted explanation of this unusual behaviour in the propagation of the 
primaries - prior to their impact with the earth atmosphere - exists yet. A large 
experimental effort [3, 4] in the next several years will hopefully clarify several of 
the issues related to this behaviour. 

While the ankle is mentioned in the debate regarding the possible existence 
of the so called Greisen, Zatsepin and Kuzmin (GZK) cutoff [5], due to the 
interaction of the primaries with the cosmic background radiation, the proposed 
resolutions of this puzzle are several, ranging from a resonant Z-burst mechanism 
[6] to string relics and other exotic particle decays [7, 8, 9, 10]. The existence of 
data beyond the cutoff has also been critically discussed [11]. 

Given the large energy involved in the first stage of the formation of the air 
showers, the study of the properties of the cascade should be sensitive to any new 
physics between the electrowcak scale and the original collision scale. Especially 
in the highest energy region of the spectrum, the energy available in the inter- 
action of the primaries with the atmospheric nuclei is far above any conceivable 
energy scale attainable at future ground-based accelerators. Therefore, the pos- 
sibility of detecting supersymmetry, for instance, in cosmic ray showers has also 
been contemplated [12]. Thus, it is not surprising, that most of the attempts to 
explain these features of the cosmic ray spectrum typically assume some form of 
new physics at those energies. 

With the advent of theories with a low fundamental scale of gravity [13] 
and large compact or non-compact extra dimensions, the possibility of copiously 
producing mini black holes (based on Thome's hoop conjecture [14]) in collisions 
involving hadronic factorization scales above 1 TeV has received considerable 
attention [15] [16] and these ideas, naturally, have found their way also in the 
literature of high energy cosmic rays [17, 18] and astrophysics [19]. For instance, it 
was recently suggested that the long known Centauro events might be understood 
as evaporating mini black holes, produced by the collision of a very energetic 
primary (maybe a neutrino) with a nucleon (quark) in the atmosphere [20]. Other 
proposals [21] also either involve new forms of matter (for example strangelets) 
or speculate about major changes in the strong interaction dynamics [22]. 

While estimates for the frequencies of these types of processes both in cosmic 
rays [23, 17, 20] and at colliders [16] have been presented, detailed studies of 
the multiplicities of the particles collected at the detectors, generated by the 
extensive atmospheric air showers following the first impact of the primary rays. 
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are far from covering all the main features of the cascade [26] . These studies will 
be useful in order to eventually disentangle new physics starting from an analysis 
of the geometry of the shower, of the multiplicity distributions of its main sub- 
components [27] and of its directionality from deep space. For instance, the study 
of the location of the maxima of the showers at positions which can be detected 
by fluorescence mirrors [28], generated as they go across the atmosphere, and 
their variations as a function of the parameters of the underlying physical theory, 
may help in this effort [23]; other obscrvables which also contain potential new 
information are the multiplicities of the various particle sub- components and the 
opening of the showers as they are detected on the ground [27]. We will focus on 
this last type of observables. 

To summarize: in the context of the TeV scale gravity with large extra di- 
mensions it is reasonable to assume that mini black holes, black holes with mass 
of a few TeV, can form at the first impact of ultra high energy primary cosmic 
rays with nucleons in the atmosphere. The black hole will evaporate into all 
types of particles of the Standard Model and gravity. The initial partons will 
hadronize and all resulting particles as they propagate in the atmosphere will 
develop into a shower (s), which eventually will reach the detectors. The nature 
and basic characteristics of these showers is the question that is the main subject 
of the present work. What is the signature on the detector of the showers arising 
from the decay of such mini black holes and how it compares with a normal (not 
black hole mediated) cosmic ray event, due, for instance, to a primary proton 
with the same energy colliding with an atmospheric nucleon (the "benchmark" 
event used here). The comparison will be based on appropriate observables of 
the type mentioned above. 

Our incomplete control of the quantum gravity /string theory effects, of the 
physics of low energy non-perturbative QCD and of the nature of the quark-gluon 
plasma phase in QCD, makes a fully general analysis of the above phenomena 
impossible at this stage. To proceed, we made the following simplifying assump- 
tions and approximations. (1) The brane tension was assumed much smaller than 
the fundamental gravity scale, so it does not modify the flat background metric. 
It is not clear at this point how severe this assumption is, since it is related to the 
"cosmological constant problem" and to the concrete realization of the Brane- 
World scenario. (2) The black hole was assumed to evaporate instantly, leading 
to initial "partons" , whose number and distributions are obtained semiclassically. 
No virtual holes were discussed and no back reaction was taken into account. (3) 
The initial decay products were assumed to fly away and hadronize, with no 
intermediate formation of a quark-gluon plasma or of a disoriented chiral con- 
densate (DCC). (4) We used standard simulation programs for the investigation 
of the extensive air showers produced in the cases of interest. To this purpose, we 
have decided to use the Monte Carlo program CORSIKA [29] with the hadronic 
interaction implemented in SIBYLL [30] in order to perform this comparison, 
selecting a benchmark process which can be realistically simulated by this Monte 
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Carlo, though other hadronization models are also available [31]. Finally, (5) 
a comment is in order about our selection of benchmark process and choice of 
interesting events. In contrast to the case of a hadronic primary, the mini black 
hole production cross section due to the collision of a > 10'^ TeV neutrino with a 
parton is of the order of the weak interaction neutrino-parton cross section [20] . 
It would, thus, be interesting to compare the atmospheric showers of a normal 
neutrino-induced cosmic ray event to one with a black hole intermediate state. 
Unfortunately, at present neutrinos are not available as primaries in CORSIKA, 
a fact which sets a limitation on our benchmark study. However, it has to be 
mentioned that neutrino scattering off protons is not treated coherently at very 
high energy, since effects of parton saturation have not yet been implemented in 
the existing codes [27]. As shown in [32] these effects tend to lower the cross 
section in the neutrino case. For a proton-proton impact, the distribution of mo- 
menta among the partons and the presence of a lower factorization scale should 
render this effect less pronounced. For these reasons we have selected as bench- 
mark process a proton-to-air collision at the same depth {Xq) and with the same 
energy as the corresponding "signal event" . In order to reduce the large statisti- 
cal fluctuations in the formation of the extensive air showers after the collisions, 
we have chosen at a first stage, in the bulk of our work, to simulate collisions 
taking place in the lower part of the atmosphere, up to 1 km above the detector, 
in order to see whether any deviation from a standard scattering scenario can be 
identified. Another motivation for the analysis of such deeply penetrating events 
is their relevance in the study of the possibility to interpret the Centauro events 
as evaporating mini black holes [20]. A second group of simulations have been 
performed at a higher altitude, for comparison. 

The present paper consists of seven sections, of which this Introduction is the 
first. In Section 2 we briefly describe the D-brane world scenario, in order to make 
clear the fundamental theoretical assumptions in our study. A brief review of the 
properties of black holes and black hole evaporation is offered here, together with 
all basic semiclassical formulas used in the analysis, with the dependence on the 
large extra dimensions shown explicitly. In Section 3 a detailed phenomenological 
description of the modeling of the decay of the black hole is presented, which is 
complementary to the previous literature and provides an independent character- 
ization of the structure of the decay. Incidentally, a Monte Carlo code for black 
hole decay has also been presented recently [33]. We recall that this description 
-as done in all the previous works on the subject- is limited to the Schwarzschild 
phase of the lifetime of the mini black hole. The modeling of the radiation emis- 
sion from the black hole - as obtained in the semiclassical picture - (see [34] for 
an overview) is performed here independently, using semi-analytical methods, 
and has been included in the computer code that we have written and used, and 
which is interfaced with CORSIKA. Recent computations of the greybody factors 
for bulk/brane emissions [34], which match well with the analytical approach of 
[35] valid in the low energy limit of particle emission by the black hole, have also 
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been taken into account. Section 4 contains our modeling of the hadronization 
process. The hadronization of the partons emitted by the black hole is treated 
analytically in the black hole rest frame, by solving the evolution equations for 
the parton fragmentation functions, making use of a special algorithm [36] and of 
a specific set of initial conditions for these functions [37]. After a brief discussion 
in Section 5 of the transformation of the kinematics of the black hole decay event 
from the black hole frame to the laboratory frame, we proceed in Section 6 with a 
Monte Carlo simulation of the extensive air showers of the particles produced by 
taking these particles as primaries. The simulations are quite intensive and have 
been performed on a small computer cluster. As we have already mentioned, in 
this work we focus on the multiplicities, on the lateral distributions of the events 
and on the ratio of electromagnetic to hadronic energies and multiplicities and 
scan the entire ultra high energy part of the cosmic ray spectrum. Our results are 
summarized in a series of plots and are commented upon in the final discussion 
Section 7. 

2 TeV Scale Gravity, Large Extra Dimensions 
and Mini Black Holes 

The theoretical framework of the present study is the D-brane world scenario 
[13]. The World, in this scenario, is 10 dimensional, but all the Standard Model 
matter and forces are confined on a 4 + dimensional hypcrsurfacc (the D^^^l' 
brane). Only gravity with a characteristic scale can propagate in the bulk. 
The TiL longitudinal dimensions are constrained experimentally to be smaller than 
0{TeV~^). However, for our purposes these dimensions may be neglected, since 
the Kaluza-Klein excitations related to these have masses at least of 0{TeV~^), 
too large to affect our discussion below. Consistency with the observed Newton's 
law, on the other hand, leads to the relation Mpi = M^+^V^, between Mpi ~ 10^^ 
GeV, the fundamental gravity scale M* and the volume Vn of the n = 6 — ul 
dimensional compact or non-compact transverse space. A natural choice for M*, 
dictated a priori by the "gauge hierarchy" puzzle, is M* = 0{Mw) — 0{1 TeV), 
while the simplest choice for the transverse space is an n-dimensional torus with 
all radii equal to R. Thus, one obtains a condition between the number n and the 
size R{n) of the transverse dimensions. Notice that under the above assumptions 
and for all values of n, it! is much larger than 10~^^cm, the length scale at which 
one traditionally expects possible deviations from the 3-dimensional gravity force, 
and the corresponding dimensions are termed "large extra dimensions" (LED). 
For n = 2 one obtains R{n = 2) of the order of a fraction of a mm. At distances 
much smaller than R one should observe 3 + n-dimensional Newton's law, for 
instance, as in torsion balance experiments [38]. Current bounds on the size of 
these large extra dimensions and on M* come from various arguments, mostly of 
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astrophysical (for instance M* > 1500 TeV for n = 2) or cosmological (M* > 1.5 
TeV for n = 4) origin [34] . A larger number of LED (n) translates into a reduced 
lower bound on Af*. It should be pointed out, that in general it is possible, even 
if "unatural" , that the transverse space has a few dimensions large and the others 
small. Here we shall assume a value of M* of order ITeV, neglect the small extra 
dimensions and treat the number of LED (n) as a free parameter. 

The implications of the existence of LED are quite direct in the case of black 
hole physics. The black hole is effectively 4-dimensional if its horizon (r//) is 
larger than the size of the extra dimensions. In the opposite case {th <^ R, 
or equivalently for black hole masses Mbh ^ lO-'^^kg for n = 6 [20]) it is 4 + n 
dimensional, it spreads over the full space and its properties are those of a genuine 
higher dimensional hole. According to some estimates, over which however there 
is no universal consensus [39], black holes should be produced copiously [15] [40] 
in particle collisions, whenever the center of mass energy available in the collision 
is considerably larger than the effective scale {^/s » M*). With ~1 
TeV, one may contemplate the possibility of producing black holes with masses 
of order a few TeV. 

Their characteristic temperature Tfj is inversely proportional to the radius 
Th of the horizon, or roughly of order M* and evaporate by emitting particles, 
whose mass is smaller than Th- The radiation emitted depends both on the 
spin of the emitted particle, on the dimension of the ambient space and on the 
amount of back-scattering outside the horizon, contributions which are commonly 
included in the so called "greybody factor", which are particularly relevant in the 
characterization of the spectrum at lower and at intermediate energy. A main 
feature of the decaying mini black hole is its large partonic multiplicity, with a 
structure of the event which is approximately spheroidal in the black hole rest 
frame. 

Once produced, these mini black holes evaporate almost instantly. The phe- 
nomenological study of 4-dimensional black holes of large mass and, in particular, 
of their Hawking radiation [41] [42] [43] , as well as the study of the scattering of 
states of various spins (s = 0,1/2,1) on a black hole background, all performed 
in the semiclassical approximation, have a long history. For rotating black holes 
one identifies four phases characterizing its decay, which are (1) the balding phase 
(during which the hole gets rid of its hair); (2) the spin-down phase (during 
which the hole slows down its rotational motion); (3) the Schwarzschild phase 
(the usual semiclassically approximated evaporation phase) and, finally, (4) the 
Planck phase (the final explosive part of the evaporation process, with important 
quantum gravitational contributions). Undoubtedly, the best understood among 
these phases is the Schwarzschild phase, which is characterized by the emission 
of a (black body) energy spectrum which is approximately thermal, with a super- 
imposed energy-dependent modulation, especially at larger values of the energy. 
The modulation is a function of the spin and is calculable analytically only at 
small energies. Extensions of these results to 4-|-n dimensions are now available. 
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especially in the Schwarzschild phase, where no rotation and no charge parameter 
characterize the background black hole solutions. Partial results exist for the spin 
down phase, where the behaviour of the greybody factors have been studied (at 
least for 1 additional extra dimension) both analytically and numerically. The 
Planck phase, not so relevant for a hole of large mass (say of the mass of the sun 
(M ~ 2 X 10^^ gr) which emits in the nano-Kelvin region, is instead very relevant 
for the case of mini-black holes, for which the separation between the mass of the 
hole and the corresponding (effective) Planck mass gets drastically reduced 
as the temperature of the hole raises and the back-reaction of the metric has to 
be taken into account. 

In the discussion below we shall use the semiclassical formulas derived for 
large black holes in the Schwarzschild phase and naively extrapolate them to the 
mini black holes as well. This is not, we believe, a severe approximation for the 
phenomena we shall discuss. As the hole evaporates, it looses energy, its mass 
decreases, its temperature increases and the rate of evaporation becomes faster. 
Thus, the hfetime of the hole is actually shorter than the one derived ignoring the 
back reaction. As we shall see below, the naive lifetime is already many orders 
of magnitude smaller than the hadronization time. This justifies the use of the 
"sudden approximation" we are making of the decay process and explains why 
the neglect of the back reaction is not severe. 

We recall that the metric of the A + n dimensional hole in the Schwarzschild 
phase is given by [44] 



dt^- 









dr'-r'dni^^, (1) 



where n denotes the number of extra spacelike dimensions, and dfll_^_^ is the area 
of the (2 -|- n)-dimensional unit sphere which, using coordinates < </? < 27r and 
< 9i < with i — 1, n + 1 takes the form 

dfll^^ = del^^ + sin^ On+i (del + sin' ( - + sin' 62 {dOl + sin' Oi d(^') )) . (2) 



The temperature Th of the black hole is related to the size of its horizon by [44] 

r« = Jii (3) 

47r rib- 
and the formula for the horizon th can be expressed in general in terms of the 
mass of the black hole Mbh and the gravity scale [44] 

For n = and M^, = Mpi ~ 10^^ GcV it reproduces the usual formula for 
the horizon {th — 2GMbh) of a 4 dimensional black hole. For n > the 
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relation between rn and Mbh becomes nonlinear and the presence of M* in the 
denominator of Eq. (4) in place of Mpi increases the horizon size for a given Mbh- 
For Mbh/M^ ~ 5 and = 1 TeV the size of the horizon is around 10~^ fm and 
decreases with increasing n. 

In the Schwarzschild/spin-down phase, the number of particles emitted per 
unit time by the black hole as a function of energy is expressed in terms of the 
absorption/emission cross sections aj^}^{uj) (or equivalently of the grcybody factors 
r{uj)), which, apart from n, depend on the spin (s) of the emitted particle, the 
angular momentum (j) of the partial wave and the corresponding energy (a;), 

dtdcu y 27r2 ex.^{u/TH)±l ^' 

Multiplying the rate of emitted particles per energy interval dN^^\u)/dtdw by 
the particle energy uj one obtains for the power emission density 

(^)/',.,^ ,.,3 



dtdu ^ 27r2 exp(a;/Tj^)±l'^'^ 

where the sum is over all Standard Model particles and the +(— ) in the denom- 
inator correspond to fermions (bosons), respectively, crj*^ are the cross sections 
for the various partial waves and depend on the spin s of each particle. We recall, 
that in the geometric optics approximation a black hole acts as a perfect absorber 
of slightly larger radius than th [45], which can be identified as the critical 
radius for null geodesies 



^ + 3^1/(n+l) /^ + 3 

rn. (7) 



V n - 

The optical cross section is then defined in function of (or equivalently 
th via Eq. (7)), such that A^, the effective surface area of the black hole hole 
projected over a k-dimensional sub-manifold becomes [46]: 



and 



2n 2 



is the volume of a k-sphere. 

It is convenient to rewrite the greybody factors as a dimensionless constant 
Ts = o-g/A^ normalized to the effective area of the horizon A4, obtained from (8) 
setting k — A and d — A + n 



^4 = 47r(^) ^r^, (10) 
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and replacing the particle cross section a in terms of a thermal averaged graybody 
factor r,(ri/2 = 2/3, Fi = 1/4, Fq — 1, i denoting the spin or species [47]). 
Eqs. (5) integrated over the frequency give (for particle i) 

^ = a{n,rH)n (H) 

with 

th) = ^ r, F(3) C(3) Q A4 r|, (12) 

where q is the number of degrees of freedom of particle i and fi is defined by the 
integral (sj is the spin) 

i -/-r(3)C(3)r3 (13) 

from which fi = 1 {fi = 3/4) for bosons (fermions). These numbers depend on 
the dimension of the brane, which in our case is 3. T{x) and ({x) are the Gamma 
and the Riemann function respectively. Since depends on the temperature 
(via Tif ) , after some manipulations one obtains 

^4^1 = -(^) (n + 3)(n + l)r^ (14) 

and 

Summing over all the particles i we obtain the compact expression 

^ = ^(E/^r,c,) F(3)c(3)r^ (16) 

with 

F,(n + l)(n + 3)(--+=^)/("+i) 
The emission rates are given by 

iV..4x3.7xlO^. '"^^'~n»^^ (-gL).-.. (18) 
21 (n + 3)<»+-'>/(»+i)(n + 1) / r„ \ _i 



^.«4x 3.7x10-^ '-^^ (19) 
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for particles with s = 0, 1/2, 1, respectively. Furthermore, integration of Eq. (6) 
gives for the black hole mass evolution 

dM _ dE 4 

- ^ (e r(4)c(4) T|, 

(21) 



with 



/5 = ^E(c^r,/',) A4r(4)c(4) (22) 



where now f[ = 1 (7/8) for bosons (fermions). Taking the ratio of the two 
equations (21) and (16) we obtain 

dN /a\ 1 



dM \(5) Th 

n + 1 \Mj 



pMMJ_fil^'™', (23) 



where we have defined 



and 



E.Q/,r,r(3)c(3) 



(25) 



E.c,Ar,r(4)C(4)- 

This formula does not include corrections from emission in the bulk. 

In the "sudden approximation" in which the black hole decays at its origi 



nal formation temperature one easily finds = {^^^^^j, where E is the energy 

spectrum of the decay products, and using Boltzmann statistics N one 
obtains the expression [15] 

N=^(^)^ein). (26) 

n+l\ M, J ^ ' ^ ^ 

This formula is approximate as are all the formulas for the multiplicities. A more 
accurate expression is obtained integrating Eq. (23) to obtain 



10 



le+06 



Iff 



10000 



. \ 



T 



O-O N(n), M=1.4 TeV 
DL(n), M=1.4TeV 
O- -O N(n), M=14 TeV 
■ - DL(n), M=14 TeV 
X-X N(n), M=141 TeV 
[>■ ■E> DL(n), M=141 TeV 




Figure 1: Multiplicities computed with Eq. (28) (N(n)) and Eq. (26) (DL(n)) for 
a varying number of extra dimensions n. DL(n) is the expression given in [15]. 



and noticing that the entropy of a black hole is given semiclassically by the 
expression 



n + lM 
n + 2Y 
Air fMBH 



one finds that 



n + 2 V M 



n+2 
n+1 



9 in), 



(28) 

which can be computed numerically for a varying n. As one can see from Fig. 1, 
the two formulas for the multiphcities are quite close, as expected, but Eq. (26) 
gives larger values for the multiplicities compared to (28) as noted by [24] . Other 
expressions for the multiplicities can be found in [25]. Since the number of ele- 
mentary states becomes quite large as we raise the black hole mass compared to 
the (fixed) gravity scale, and given the (large) statistical fiuctuations induced by 
the formation of the airshower, which reduce the dependence on the multiplicity 
formula used, we will adopt Eq. (26) in our simulations. Overall, in the masslcss 
approximation, the emission of the various species for a 3-brane is characterized 
by approximately 2% into spin zero, 85% into spin half and 13% into spin one 
particles, with similar contributions also for the power emissivities. These num- 
bers change as we vary the dimension of the brane (d) and so does the formula 
for the emissivities, since the number of brane degrees of freedom {ci{d)) has to 
be recomputed, together with the integrals on the emission spectra {fi{d)) [23]. 
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The integration of the equation for the power spectrum, in the massless ap- 
proximation, can be used to compute the total time of decay (assuming no mass 
evolution during the decay) 

which implies that at an energy of approximately 1 TeV the decay time is of the 
order of 10~^^ seconds. Therefore strong interaction effects and gravity effects 
appear to be widely separated and hadronization of the partons takes place after 
their crossing of the horizon. The black hole is assumed to decay isotropically 
(s-wave) to a set of N elementary states, selected with equal probability from all 
the possible states available in the Standard Model. We mention that in most 
of the analysis presented so far [18, 17, 23] the (semiclassical) energy loss due to 
bulk emission has not been thoroughly analyzed. We will therefore correct our 
numerical studies by keeping into account some estimates of the bulk emission. 



3 Modeling of the Black Hole Decay 

The amount of radiation emitted by the black hole in the ED is viewed, by an ob- 
server living on the brane, as missing energy compared to the energy available at 
the time when the black hole forms. From the point of view of cosmic ray physics 
missing energy channels imply reduced multiplicities in the final air shower and 
modified lateral distributions, these two features being among the main observ- 
ables of the cosmic ray event. However, since the initial energy of the original 
cosmic ray is reconstructed by a measurement of the multiplicities, an event of 
reduced multiplicity will simply be recorded as an event of lower energy. It is 
then obvious that an additional and independent reconstruction of the energy of 
the primary cosmic ray is needed in order to correctly identify the energy of these 
events. 

In our study we will compute all the observables of the induced air shower 
using both the lab frame (LF) and the black hole frame (BHF) to describe the 
impact and the formation of the intermediate black hole resonance. Also, in 
the simulations that we will perform, the observation level at which we measure 
the properties of the air showers will be selected to take properly into account 
the actual position of a hypothetical experimental detector. The target of the 
first impact of mass M is assumed to be a nucleon (or a quark) at rest in the 
atmosphere and the center of mass energy, corrected by emission loss in the bulk, 
is made promptly available for an instantaneous black hole formation and decay. 
We will also assume that the energy E\ of the incoming primary varies over all 
the highest part of the cosmic ray spectrum, from 10^^ eV up to lO^'^ eV. 

We denote by /3 the speed of the black hole in the lab frame. In our notations, 
E* is the typical energy of each elementary state in the decay products (parton. 
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lepton) in the BH frame and P* is its corresponding momentum. 

We will assume that a black hole decays "democratically" into all the possible 
partonic states, proportionally to the number of Standard Model states which are 
available to it at a given energy. 

The energy per partonic channel will be appropriately weighted and we will 
assume that each parton (/) will decay into a final state hadron h (carrying a 
fraction x of the original momentum), with a probability distribution given by 
the corresponding fragmentation function D^{x, Q"^), which is evolved from some 
low energy input scale Qo up to the relevant scale characterizing the decay. This 
is given by the available energy per fundamental state, equally distributed among 
all the states. 

The quantification of the injection spectrum involves a computation of the 
relevant probabilities for the formation of all the possible hadronic/leptonic states 
prior to the simulation of the air shower. Let's briefly elaborate on this. 

To move from the parton level to hadron level, we let D'^{x, Q^), D^{x, Q'^), 
and Dg{x, Q"^) be the fragmentation functions of Np quarks 5, antiquaries and 
of the gluon g, respectively, into some hadron h with momentum fraction x at 
the scale Q. From the fragmentation functions we obtain, for each hadron h, 
the mean multiplicity of the corresponding s-wave and the corresponding average 
energy and momentum. Specifically we obtain 

<Dn> = E r dzD'}{z,Q^) (30) 
for the probability of producing a hadron h, and 

E*H = ET zdzD'}{z,Q') (31) 

for the average energy of the same hadron. We recall that Zmin is the minimal 
fraction of energy a hadron (/i), of mass rrih, can take at a scale Q, and can be 
defined as z^in — mh/{Q/2). In practical applications one can take the nominal 
value Zmin = 0.05 for every hadron, without affecting much the mean multiplici- 
ties and the related probabilities. This implies that 

< > + E < > + < > + < >= Prh (32) 
/ 

together with the condition J2h = 1, where the sum runs over all the types of 
hadrons allowed by the fragmentation. In all the equations above, the fragmen- 
tation takes place at the typical scale Q = E/N, scale at which the moments are 
computed numerically. For the identification of the probabilities it is convenient 
to organize the 123 fundamental states of the Standard Model into a set of flavour 
states (g/), with / running over all the flavours except for the top quark, where 
in (g/) we lump antiquark states and color states, plus some additional states. 
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The weight of the (qf) setisp/ = 2x2x3/123, where the factors 2 and 3 refer 
to spin, quark-antiquark degeneracy and color. It is worth to recall that quark 
and antiquark states of the same flavour have equal fragmentation functions in 
all the hadrons, and this justifies the q/q degeneracy of the set. The additional 
states are the gluon (g) with a weight Pg ~ 2 x 8/123, the photon (7), with a 
weight Pj — 2/123 and the remaining states (r) in which we lump all the states 
which have been unaccounted for, whose probabilities Pr are computed by differ- 
ence. These include the top and the antitop (12/123), the W's and Z (9/123) and 
the leptons (24/123). The fragmentation functions into hadrons, corresponding 
to these states, < > are computed by difference from the remaining ones 
< Dg >, < Df > and < >, which are known at any scale Q from the litera- 
ture. Beside the favour index f = u,d, c, s, b, introduced above, we introduce a 
second index i running over the (r) states, the photon and the gluons {i = g, 7, r). 

The probability of generating a specific sequence of states in the course of 
the evaporation of the black hole is then given by a multinomial distribution of 
the form 

finf,n,,Pf,p,) = n .Tn^l IIPyilP^"' (33) 

which describes a typical multi-poissonian process with N trials. Notice that, to 
ensure proper normalization, we need to require that 

JJ^j! = ng\n^\nr\ 

i 

= ng\n^\{N -Y^rif -ng-n^)\ (34) 
/ 

The computation of the cumulative probabilities to produce any number of 
hadrons of type h by the decay of the black hole are obtained from the multino- 
mial distribution multiplied by the fragmentation probabilities of each elementary 
state into {h) and summing over all the possible sequences 

ftcum h - E n Jn „, n (pf < D) >)"' n (p. < A" >)"• . (35) 

A possible way to compute Pr^^^j^ when iV is large is to multiply the multi- 
nomial distribution by a suppression factor Exp[— A(5]]j + J^j^f — N)], with 
A a very large number, and interpret this factor as a Boltzmann weight, as in 
standard Monte Carlo computations of the partition function for a statistical 
system. Simulations can be easily done by a Metropolis algorithm and the con- 
figurations of integers selected are those for which the normalization condition 
N = J2i f^i + J2 f ''T' f is satisfied. In our case, since we are interested only in the 
mean number of hadrons produced in the decay and in their thermal spectrum, 
the computation simplifies if we average over all the relevant configurations. 
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4 Fragmentation and the Photon Component 



The evolution with of the fragmentation functions is conveniently formulated 
in terms of the linear combinations 



Np 
i=l 



(36) 



(37) 



(38) 



as for these the gluon decouples from the non-singlet (+) and the asymmetric 
(-) combinations, leaving only the singlet and the gluon fragmentation functions 
coupled; 



dQ^ 



(39) 
(40) 



PE(«s(Q'))<8)i?l(g')] {x) 



+2Np 



P,^G{as{Q'))®D^a{Q')]{^): (41) 



2Nf 



Pg^, {as{Q^)) ® D^{Q^)\ {x) 
P,^, {a^m) Dj(Q^)] {x) . 



+ 



The kernels that appear in the equations above are defined by 

P(+) {x,as{Q'')) = Pl^^{x,as{Q^))+Pl^-^{x,as{Q')) , 

Pi:{x,as{Q^)) = P(+) {x,as{Q^))+2NFPgt,{x,as{Q^)) , 

P(_) {x, as{Q')) = {x, as{Q')) - {x, a,{Q^)) , 



(42) 

(43) 
(44) 
(45) 



with as{Q^) being the QCD coupling constant. In the perturbative expansion of 
the splitting functions, 



P{x,a,{Q')) = 4^p(°)(x) + ( ^Vp«(x) + ^ ^"^^^ 



27r 



27r 



27r 



■ (46) 
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The timelike kernels that we use are given by 

"3 



G^G 



X 



Cf 



5{l-x) + 2(^-^ 
0, 



X 



— 1 — X 



i + {i-xy 



2NfTr 



X 



x^ + (l-xf 



6 ^ 



-NfTr]S{1-x) 



+2Nc 



1 \ 1 
+--2 + X 

1 — X/ + X 



X 



(47) 
(48) 
(49) 
(50) 

(51) 



The formal solution of the equations is given by 

(52) 

where Qq is the starting scale of the initial conditions, given by D[x,Qq). At 
leadiong order in as, we solve this equation using a special ansatz 

'-^■^)^" (53) 



C}(x,Q') = i:^log 



ni 



MQl) 



and generating recurrence relations at the n + l — th order for the An+i coefficients 
in terms of the [36] . It is easy to see that this corresponds to the numerical 
implementation of the formal solution 

D){x, = Exp {tP®) D){x, Ql) (54) 

with t = {oisiQ"^) / OisiQo)) , where the exponential is a formal expression for an 
infinite iteration of convolution products. We show in Figures 2-4 results for 
some of the fragmentation functions into pions, kaons and protons computed for 
a typical parton scale of 200 GeV. 

The photon contributions to the decay of the black hole is treated separately. 
The evolution equation for the fragmentation functions of photons and parton 
fragmentation into photon D^(a;, Q^), D'^{x,Q^) satisy at leading order in ctem 
(the QED fine structure constant) and ag (the QCD coupling), the evolution 
equations [49] 



dlng2 

which can be integrated with the initial conditions D'^{x, Q"^) — 5{1 — x), and 



(55) 



dlng2 



a 



^—P,^,{x)(^D;{x,Q') 



(56) 
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Figure 4: Fragmentation functions into p/p at 200 GeV. 



which can also be integrated with the resuh [50] 



(57) 



In Eq. (57) the second term is termed the hadronic boundary conditions, which 
come from an experimental fit, while the first term is the pointlike contributions, 
which can be obtained perturbatively. In [51] a leading order fit was given 



a 



Dl{x, Qi) = — [- P,^,{x) ln(l - xY - 13.26 



(58) 



at the starting scale Qq = 0.14 GeV. The kernels in these cases arc given by simple 
modifications of the ordinary QCD kernels, for instance Pq^j{x) — Cq/CpPq^gix) 
[49]. 

For the fragmentation function of quarks to photons with virtuality My, the 
perturbative result is given in [52] 



a 
2^ 



l + il-x)\ xQ^ / P2 

In -^t:: X \ 1 — 



X 



M2 



(59) 



where is the virtuality of the photon. The gluon to photon transitions are 
neglected, since Pg-*^ vanishes in leading order. 

We recall that each elementary state emitted is characterized by an average 
energy given by (e) = Mbh/ {N). 

The leptonic component e^, //^, produced by the decay is left unaltered and 
provides an input for the air shower simulator as soon as these particles cross the 
horizon. The r^s are left to decay into their main channels, while the hadroniza- 
tion of the u, d, s, c quarks and the gluons is treated with our code, that evolves 
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11 


0.451 


0.226 


0.048 


0.174 


0.067 


0.034 








U.UOt: 








d 


0.451 


0.226 


0.174 


0.048 


0.034 


0.067 






231 


252 


084 


03^ 


070 


s 


0.391 


0.195 


0.068 


0.068 


0.139 


0.139 




0.295 


0.117 


0.081 


0.081 


0.108 


0.108 


c 


0.329 


0.165 


0.167 


0.167 


0.085 


0.085 




0.309 


0.155 


0.194 


0.194 


0.071 


0.071 


b 


0.438 


0.219 


0.129 


0.129 


0.042 


0.042 




0.324 


0.162 


0.115 


0.115 


0.041 


0.041 


9 


0.303 


0.152 


0.253 


0.253 


0.020 


0.020 




0.807 


0.404 


0.317 


0.317 


0.034 


0.034 



Table 1: Initial conditions. For each couple of parton and hadron, the upper 
number in the box is the probability for the parton / to hadronize into the 
hadron h, (/J^.^ D){z, Q)dz) / Eh' ^/'(^, Q)dz), while the lower number is 

the average energy fraction of h, J^,^-^ zDj{z,Q)dz. In this table the energy of 
the parton is Q = 1.414 GeV for u, d, s and g, Q = 2mc = 2.9968 GeV for c and 
Q = 2mh = 9.46036 GeV for the b quark, generated via the set of ref. [37] 

the fragmentation functions to the energy scale (s). The top (t) quark is treated 
consistently with all its fundamental decays included; hadronization of the b 
quark is treated with suitable fragmentation functions and also involves a suit- 
able evolution. As we vary Mbh and wc scan over the spectrum of the incoming 
cosmic rays the procedure is repeated and rendered automatic by combining in a 
single algorithm all the intermediate steps. Tables 1 and 2 contain the results of a 
renormalization group analysis of the fragmentation functions for all the partons 
(except the top quark), where we show both the intial conditions at the input 
scale, whose lowest value is Q = 1.414 GeV, and the results of the evolution, at 
a final scale of Q = 200 GeV, the initial set being taken from ref. [37]. 

This concludes the computation of the probabilities for each hadron/lepton 
present in the decay products of the mini black hole. It is reasonable to assume 
that these particles will be produced spherically, since higher angular momenta 
are suppressed by the corresponding centrifugal barrier. However, the analysis 
of the shower profile has to be performed in the lab frame. This requires the 
transformation of the initial configurations above to the laboratory frame, which 
is exactly what is discussed next. 
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0.446 


0.223 


0.079 


0.166 


0.053 


0.033 
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077 


1 78 


047 


097 
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0.446 


0.223 


0.166 


0.079 


0.033 


0.053 






193 


1 78 


077 


027 


047 


s 


0.425 


0.213 


0.093 


0.093 


0.088 


0.088 




0.295 


0.117 


0.077 


0.077 


0.070 


0.070 


c 


0.371 


0.185 


0.158 


0.158 


0.064 


0.064 




0.295 


0.147 


0.150 


0.150 


0.051 


0.051 


b 


0.431 


0.216 


0.132 


0.132 


0.045 


0.045 




0.292 


0.146 


0.101 


0.101 


0.036 


0.036 


9 


0.428 


0.214 


0.135 


0.135 


0.044 


0.044 




0.577 


0.289 


0.175 


0.175 


0.057 


0.057 



Table 2: For each couple of parton/hadron, the first number is 
the probability of fragmentation of the parton / into the hadron 
(■tin ^/(^' I ^h' {izmin ^/'(^' Q)^^) , whllc thc sccoud is the average en- 

ergy fraction of h, Jl^.^ zD^{z, Q)dz. The energy of the parton is Q = 200 GeV. 



5 Sphericity and Boost 

The transformation from the black hole frame (BHF) to the laboratory frame 
(LP) is performed by a Lorentz boost with speed /3, the speed of the black hole in 
the LF. Assuming that the black hole is produced in the collision of a primary of 
energy Ei in the LF and negligible mass compared to Ei^ with a parton of mass 
M in the atmosphere, one obtains j3 = Ei/{Ei + M). A spherical distribution of a 
particular particle of mass m among the decay products in the BHF is transformed 
to an elliptical one, whose detailed form is conveniently parametrized by 



P 1- 



M 



where f3* — P*/E* is the speed of this particle in the BHF, the ratio of its BHF 

momentum and energy. Figure 5 depicts the relevant kinematics. A particle 
emitted in the direction 6* in the BHF, is seen in the direction 6 in the LF, with 

tang^yr^ ^f^* (61) 
* g* + cos 0* 

For values of 5^* > 1 the shape of the 1-particle distribution in the LF is charac- 
terized by a maximum angle of emission 




tan Omax 1 = 1/ —5 7 , (62) 
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Figure 5: An s-wave distribution in the BHF is transformed in the detector frame 
to an eUiptical one, whose detailed shape depends on the value of g* 



which for ^f* = 1 is equal to 90°. Only for g* < 1 there is backward emission 
in the LF. As a relevant example, let us consider the case of a hadron of mass 
m — 1 GeV and energy E* — 100 GeV emitted by a black hole, formed by an 
initial primary of energy Ei — 1000 TeV, which hit a quark of mass M ~ lOMeV 
to form a black hole. The gives g* = 1.00005 and the corresponding maximum 
angle in the LF is 9max — ^^'^^max = 1-4 x 10~^, giving an angular opening of the 
decay products of about 2 degrees. 

As shown in Figure 5, for g* > 1, which is relevant for our purposes, the 
mapping from 9* to 9 is not one to one. In a given direction 9 in the LF, one 
receives particles emitted in two different directions 9*^ in the BHF. They satisfy 
[48] 

c^cos^ ~ [l^) (±cos^x/Z) ^ ' 

where 

K = l + -'/^{l- g*^)tan^9, (64) 
and with the momenta and energies E^ of the two branches given by 

, P*cos9{g* ±^/K) 



and 

m[j* ±vcos9 (v*^7*^ - t;^7^ sin^ 9^ 



7(1 — cos^ 9) 

1/2 



(65) 



^ 7(l-^;2cos2^) ^^^^ 
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respectively. In the above formulas 7""^ = \/\ — v{v*^ is the speed of the 
hadron in the LF (BHF), and (7*)"^ = Tn/El^ — \/l — v*'^. For massless final 
state particles, in particular, these relations become 

p* 

P = E = — (67) 

7(1 — cost/) 

and reduce to the familiar Doppler formula when ^ = 0. 

The probability distribution W^{cos6* ,(f)) of a hadron (h) as a function of 
the direction Q = {cos6*,(f)) in the BHF, assumed spherically symmetric and 
normalized to the total probability Pvh of detecting this hadron among the decay 
products with elementary states, is 

W^(cose*) = ^. (68) 



The corresponding one in the LF is 

dcos9 



W,{cose)=Y.-^^W*{cose*^). (69) 



± 

In the special case g* — 1, the probability distribution simplifies to 

COS 

Wnicos e) = 2Pr, -— , (70) 

7"'(1 — p^cosb'"')"' 

peaked in the forward direction, symmetric around the maximum value, obtained 
for = and equal to 27^, while the momentum distribution is 

T^(a\ /^7* cos^ 

(71) 

As we have already mentioned, the structure of the partonic event (and, simi- 
larly, of the hadronic event after fragmentation) is characterized by the formation 
of an elliptical distribution of partons, strongly boosted toward the detector along 
the vertical direction. Each uniform (s-wave) distribution is strongly elongated 
along the arrival direction (due to the large speed of the black hole along this 
direction) and is characterized by two sub-components (PF^), identified by a ± 
superscript. Their sum is the total probability distribution given in (63). The 
momentum component is largely dominant and strongly peaked around the 
vertical direction with rather small opening angles and this behaviour can be 
analyzed numerically with its N dependence. In the explicit identification of the 
two independent distributions ± in terms of the opening angle 0, as measured in 
the LF, we use the relations 

^y±(^) = lly±(cosr)|^|^^sin^, (72) 

z a cos y 
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where we have introduced a factor of 1/2 for a correct normahzation of the new 
distribution in the 6 variable. In Figures 6 and 7 we show the structure of these 
distributions in the LF. Both are characterized by a very small opening angle (6) 
with respect to the azimuthal direction of the incoming cosmic ray, being the 
dominant one. Two similar plots (Figures 8 and 9) illustrate the two components 
as functions of the same angle. 

One can easily check that we can now integrate symmetrically on both dis- 
tributions to obtain the correct normalization (to Pr/j) for a given hadron (or 
parton) 

W^{9)d9= / W+{9)d9+ / W-{9)d9+ / W-{9)d9^Prh (73) 
or, equivalently, using cos ^ as a distribution variable 

d rns 9* 

W^(cos9) = W^(cos9*) I ^^f^ I sin^ (74) 

dcosO 

with 

W"^ {cos 9)d cos 9 = r W+{cos9*)dcos9*+ r W- {cos 9*)d cos 9* ^PTh 

JcOsOmax JcOsOl* J — 1 

(75) 

and 9*1 obtained from (62). 
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Figure 6: Plot of the branch of the probabihty distribution for an incoming 
energy of the cosmic ray Ei = 10^ GeV and for various values of N of the 
elementary partonic states emitted during the decay of the black hole 
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Figure 7: Same plot as above but for the W branch of the distribution 
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Figure 8: Plot of the branch of the momentum distribution in eV for an 
incoming energy of the cosmic ray Ei = 10^ GeV and for various values of of 
the elementary partonic states emitted during the decay of the black hole 
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Figure 9: Same plot as above but for the P branch of the momentum distribution 
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n 


/. 





0.70711 


1 


0.66533 


2 


0.63894 


3 


0.62057 


4 


0.60696 



Table 3: Fraction of Ecm that is bound into the black hole as a function of 
the number n of extra-dimension in head-on collisions. 

6 Air Shower Simulations 

The simulation of the events is performed at the last stage, using an air shower 
simulator. We have used CORSIKA [29] with appropriate initial conditions on 
the spectrum of the incoming particles in order to generate the full event mea- 
sured at detector level. In most of the simulations we have assumed that the first 
impact takes place in the lower part of the atmosphere, not far from the level 
of the detector, at a varying altitude. The reason is that one of our interests 
is the investigation of the possibility that the Centauro events may be related 
to evaporating mini black holes, formed by the collision of weakly interacting 
particles (e.g. neutrinos) which penetrate the atmosphere. Of course, we have 
simulated events happening at higher altitudes as well. We have performed two 
separate sets of simulations, the first set being benchmark events with an "equiv- 
alent" proton replacing the neutrino-nucleon event, colliding at the same height, 
the second being the signal event, i.e. the black hole resonance. The difference 
between the first and the second set is attributed to the different components of 
the final state prior to the development of the air shower. 

We compute the average number of particles produced in the process of BH 
evaporation using the formula 

Mbh = EcMfn (76) 

where Ecm is the energy in the center of mass frame in the neutrino-nucleon 
collision and /„ is the fraction of Ecm that is bound into the black hole as a 
function of the number n of extra-dimensions. Numerical values for in head- 
on collisions are taken from Ref. [53] and reported in Table 3. 

The overall shower, defined as the superposition of the various sub-components, 
develops according to an obvious cylindrical symmetry around the vertical z-axis 
near the center. We assume in all the studies that the incoming primary under- 
goes a collision with a nucleon in the atoms of the atmosphere at zero zenith with 
respect to the plane of the detector. 

The model of the atmosphere that we have adopted consists of N2, O2 and Ar 
with the volume fractions of 78.1%, 21.0% and 0.9% [54]. The density variation 
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Layer /' 


Altitude // [km] 


(I, [g/ciir] 


b, [g/cin^] 


c, [cm] 


1 


0...4 


-186.5562 


1222.6562 


994186.38 


2 


4. ..10 


-94.919 


1144.9069 


878153.55 


3 


10... 40 


0.61289 


1305.5948 


636143.04 


4 


40 . . . 100 


0.0 


540.1778 


772170.16 


5 


> 100 


0.01128292 


1 


109 



Table 4: Parameters of the U.S. standard atmosphere. 



of the atmosphere with altitude is modeled by 5 layers. The pressure p as a 
function of the altitude h is given by 

p{h) = ai + biexp{-h/ci), i = l,...,4 (77) 

in the lower four layers and by 

p{h) ^a^-h- (78) 

in the fifth layer. 

The Oj, 6j, Q parameters, that we report in Table 4, are those of the U.S. stan- 
dard atmospheric model [29]. The boundary of the atmosphere in this model is 
defined at the height 112.8km, where the pressure vanishes. In Figure 10 we 
show a plot of the pressure {p — Xy, also called vertical depth) as a function of 
the height. 



This is defined via the integral 

/ p{h')dh' (79) 

J h 

of the atmospheric density p{h) for zero zenith angle, while the corresponding 
slant depth is given by 

roo ( 1/2 \ 

\ pNcos^+-— sin^^ (80) 
h \ 2 ity j 

for a zenith angle Q and Rt is the radius of the earth. 

To put into perspective our Monte Carlo study it is convenient to briefly 
summarize the basic features of the theory of cascades on an analytical ground. 
The theory consists of the system of transport equations [55, 56] for the numbers 
Nn{E, X) of particles of type n with energy E at height X 

dN^{E,X) Nn(E,X) 1 



dX Xn{E) CTn^PAh 



Nn{E,X), (81) 
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Figure 10: Pressure versus altitude for the U.S. standard atmosphere model. 



where Xn{E) is their interaction length, is their lifetime and 7 the Lorentz- 
factor corresponding to their given energy. In the simple case of an isothermal 
atmosphere pAir = Po exp(— /i//io) — X/ho, at a. scale height ho 



dNn{E,X) _ N^{E,X) _ 1_ 

dx - K{E) d^''"^""'^^ 



(82) 



where dn is their decay length, defined by 

1 mc^hf] 



dr, 



EcT„X 



(83) 



Particles produced at higher energies are also accounted for by an additional term 
in the cascade 



dN4E,X) 
dX 



-Nn{E,X) 



+ 



+ ^/ N^{E',X) 



XniE) dn{E)_ 

Wrnn{E',E) 
\m{E') 



(84) 



J / 7-i/\ Emn 



{E',E) 



dE' 



describing the change in the number of particles of type n due to particles of 
type m by interaction or decay, integrated over an allowed interval of energy. 
The functions Wmn{E',E) are the energy-spectra of secondary particles of type 
n in a collision of particle m with an air-molecule, while Dmn{E', E) are the cor- 
responding decay-functions. The advantage of a transport equation compared to 
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a Monte Carlo is, that it provides a rather simple analytical view of the develop- 
ment of the cascade across the atmosphere. Most common in the study of these 
equations is to use a factorized ansatz for the solution N{E,X) = A{E)B{X), 
which assumes a scaling in energy of the transition functions [55]. In our case an 
analytical treatment of the cascade corresponds to the boundary condition 



with Pr„(i?) being the probability that the black hole decays into a specific state 
n. As we have already discussed above, these decays are uncorrelated and Eq. (85) 
is replicated for all the elementary states after hadronization. The emission prob- 
abilities Pin{E) have been computed by us for a varying initial energy E = JMbh 
using renormalization group equations as described before, having corrected for 
energy loss in the bulk. The interactions in the injection spectrum of the original 
primaries at our Xq (Xq = 517 g/cm^) has been neglected since this is not im- 
plemented in CORSIKA. The showers have been performed independently and 
the results of the simulations have been statistically superimposed at the end 
with multiphcities computed at detector level {Xi — 553 g/cm^). We have kept 
the gravity scale constant at 1 TeV and varied the mass of the black hole 
according to the available center of mass energy E. As we have already discussed 
in the previous sections, as a benchmark process we have selected a proton-air 
impact at the same Xq with the boundary condition 



which occurs with probability 1 {Pvp — 1). 

We are interested both in the behaviour of the multiplicities and in the lateral 
distributions of the cascades developed at detector level. For this purpose we have 
defined the opening of the conical shower after integration over the azimuthal 
angle, as in [27], and given the symmetry of the event, we plot only the distance 
from the center as a relevant parameter of the conical shower. 

A varying number of extra dimensions n = 0, 1, . . . , 4 implies a different ratio 
for bulk-to-brane energy emission, a different average number of elementary de- 
caying states and different energy distributions among these. We have varied the 
energy Ei of the primary, thereby varying the mass of the black hole resonance, 
in the interval 10^^ — lO^'^ eV. The hadronization part of our code has been done 
by changing the PTh{E) obtained from a numerical solution of the fragmentation 
functions separately for each value of the energy shared. 

• Preliminary studies 

We start with the numerical study of the partial and total multiplicities of 
the various sub-components and of their lateral distributions in the bench- 
mark event. The results for photons and leptons are shown in Figures 11 



Nr,{E, Xo) = Pr„(£;) S{E - fMsH/iN)) 



(85) 



Np{E,Xo)^PTpS{E-fMBH) 



(86) 
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Figure 11: Benchmark event: Multiplicities of photons, e^, /x^ at an observation 
level of 4500 m as a function of the altitude of the first impact. 

to 14 as functions of the altitude of the impact and for two different obser- 
vation points at 4,500 m and 5,000 m, approximately the altitudes of the 
detectors at Pamir and Chacaltaya, respectively. These preliminary plots, 
based on actual simulations of a proton-to-air nucleus impact at 10^^ eV, 
show a steady growth of the multiplicities of the secondaries as we raise 
the point of first impact above the detector. The statistical errors in 
the simulations (80 uncorrelated events have been collected per point) are 
rather small, quite uniformly over all the altitudes of the impact, and in- 
dicate a satisfactory stability of the result. The positions of the detectors 
do not seem to have an appreciable impact on the characteristics of the 
secondaries. As for the lateral distributions we observe an increase in the 
opening of the showers with the event altitude, which is more enhanced for 
the muonic component and for the photons and less for the electrons and 
positrons. Also in this case the statistical fluctuations are rather small. On 
the basis of these results we have selected for the remaining simulations 
a first impact at 5,500 m and the observation (detector) level at 5,000 m. 
However, for comparison, we will also show later the results of a second set 
of simulations that have been performed with the first impact at 15,000 m. 

• Choice of scales and corrections 

To compare standard and black holes events, we have selected a gravity scale 
of = ITeV and varied the black hole mass, here taken to be equal to the 
available center of mass energy during the collision. Therefore, a varying 
El is directly related to a varying Mbh and we have corrected, as explained 
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Figure 12: As in Figure 11, but at an observation level of 5000 m. 
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Figure 13: Average radial opening of the shower of photons, e^, fi^ at an obser- 
vation level of 4500 m as a function of the altitude of proton's first interaction. 
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Figure 14: As in Fig. 13, but at an observation level of 5000m. 

above, for the energy loss into gravitational emission. Unfortunately, this 
can be estimated only heuristically, with bounds largely dependent on the 
impact parameter of the primary collision. A reasonable estimate may be 
of the order of 10 — 15% [35]. Corrections related to emission in the bulk 
have also been included, in the way discussed in previous sections. 

• Energy ratios: electromagnetic versus hadronic 

Not all observables are statistically insensitive to the natural fluctuations 
of the air showers. In the study of black hole versus standard (bench- 
mark) events, the study of the ratios Ai'em/A^hadron and -Ecm/-^hadron have 
been proposed as a way to distinguish between ordinary showers and other 
extra-ordinary ones. Centauro events, for instance, have been claimed to 
be characterized by a rather small ratio of electromagnetic over hadronic 
energy deposited in the detectors, contrary to normal showers, in which 
this ratio is believed to be -EcmZ-^hadron ~ 2. Instead, as one can easily 
recognize from the results presented in Figures 15 and 16, the multiplicity 
ratio takes values in two different regimes. In the "band" of values 1 — 5 
for the case of the lower first impact and 100 — 160 for the higher impact. 
The larger values of the band in this latter case are justified by the fact 
that the shower is far more developed, given the altitude of the impact, 
and therefore is characterized by an even more dominant electromagnetic 
component. The energy ratio, on the other hand can take small values, 
in agrement with the values observed in Centauros. However, notice that 
both the black hole and standard simulations show a complex pattern for 
these ratios and in addition they are characterized by large fiuctuations for 
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varying energy and number of extra dimensions. Furthermore, there does 

not seem to be a statistically significant difference in these observables be- 
tween the benchmark event and the black hole mediated ones, at least for 
event heights greater than 500 m above the detector. We conclude that (a) 
either these observables may not be suitable, especially given the limited 
statistics of the existing and future experiments, to discriminate between 
black hole mediated events versus standard ones, or (b) that the differences 
in these ratios appear in events with initial impact in the range — 500 m 
from the detector. This latter issue, of relevance to the Centauro events, is 
currently under study [57]. 

• Multiplicities 

For the multiplicities themselves the situation is much cleaner. In Fig- 
ures 17-22 we show the behaviour of the total as well as of some partial 
multiplicities in black hole mediated events and in standard events as a 
function of the energy and for a varying number of extra dimensions. 

The curves are very well fitted in a log-log plot by a linear relation of the 
form 

N = 10«(")£;i^^") (87) 

with intercepts q{n) and slopes cr{n), that increase with the number of 
extra dimensions n. We present in Figures (a) results of the simulations 
performed with a first impact taken at 5,500 m, while Figures (b) refer 
to a first coUision at 15,000 m. In Figures of type (a) the slopes of the 
benchmark events are smaller than those of the black hole events and show 
a larger intercept. This feature is common to all the sub-components of 
the air showers. A simple explanation of this fact is that at lower value of 
the impact energy, the number of states available for the decay of the black 
hole is smaller than the number of partonic degrees of freedom available in 
a proton-proton collision. We recall, as we have already discussed in the 
previous sections, that our benchmark results define in this case an upper 
bound for the total multiplicities expected in a neutrino-proton collision. 
Therefore, in a more realistic comparison, we would discover that the black 
hole and the standard results should differ more noticeably. The large 
multiphcity of the states available for the decay of the black hole dominates 
over that of a standard hadronic interaction, and this justifies the larger 
multiplicities produced at detector level. As we increase the altitude of the 
impact, in plots of type (b) we find a similar trend but the differences in the 
total and partial multiplicities are much harder to discern for black holes 
and benchmark events. In fact, for collisions starting at higher altitudes 
the showers are all fully developed and the differences between the two 
underlying events are less pronounced. 
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Figure 15: Ratio between N^m (total multiplicity of photons and e^) and N^ad 
(total multiplicity of everything else) as a function of E^. The first interaction 
is kept fixed at 5500m (517 g/cm^) (a), or at 15000m (124 g/cm?) (b), and the 
observation level is 5000 m (553 g/cw?). We show in the same plot the benchmark 
(where the primary is a proton) and mini black holes with different numbers of 
extra-dimensions n. 
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Figure 16: As in Figure 15, but this time we show the ratio between Eem (total 
energy of photons and e^) and Eh^d (total energy of everything else) as a function 
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Figure 17: Plot of the total particle multiplicity as a function of Ei. Case (a) is 
for an impact point of 5,500 m and case (b) for 15,000 m. 
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Figure 18: Plot of the multiplicity of photons as a function of (a) 5,500 m, 
(b) 15,000 m 
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Figure 19: Plot of the multiplicity of e as a function of Ei^ (a) 5,500 m, 
15,000 m 
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Figure 20: Plot of the multiplicity of as a function of E^^ (a) 5,500 m, 
15,000 m 
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Figure 21: Plot of the multiplicity of tt as a function of Ei, (a) 5,500 m, 
15,000 m. 
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Figure 22: Plot of the multiplicity of protons as a function of ^^i, (a) 5,500 m, 
(b) 15,000 m. 
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Another feature of the black hole events is that the slopes and the intercepts 
of the various plots, for a given choice of altitude of the impact, are linearly 
correlated. To illustrate this point we refer to Figures 23-24 from which 
this behaviour is immediately evident. To generate each of these figures we 
have plotted the parameters {a, q) of a corresponding plot - for the total 
or for the partial multiplicities - independently of the specific number of 
extra dimensions. The results shown in these figures clearly indicate that 
the relation between the intercept q and the slope a appearing in Eqn. (87) 
is linear and independent of n 

q^aa + P (88) 

with a and /? typical of a given setup (photons, total multiplicities, etc.) 
but insensitive to the parameter n. Therefore, black hole events are char- 
acterized by particle multiplicities on the ground of the form 

Nground = lO^'^+^^r- (89) 



• Lateral distributions 

In Figures 25-29 we illustrate the results of our study of the lateral distri- 
butions for the total inclusive shower and the various sub-components as 
a function of the incoming energy Ei. The average opening of the shower 
as measured at detector level is plotted versus energy in a log-log scale. 
Notice a growing opening of the shower as we raise the energy of the black 
hole resonance, which is more remarked for a lower number of extra di- 
mensions. In contrast, the benchmark simulation shows a small decrease 
(negative slope) with energy. The larger opening of the shower in black hole 
mediated events - compared to standard air showers - is due to the s-wave 
emission typical of a black hole decay, which is very different from an or- 
dinary collision. Contrary to the case of multiplicities, here simulations of 
type (a) and (b) show a similar trend, with very distinct features between 
standard and black hole events. Notice that in this case the difference in 
the partonic content of the two different events (benchmark versus black 
hole mediated) is less relevant, since it is the geometrical fireball emission in 
the black hole case which is responsible for the generation of larger lateral 
distributions. 

Also in this case we discover a linear relation between average radius R 
of the conical openings and energy, relation that can be fitted to a simple 
power law 

R^lOl'in)E^'^''\ (90) 
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Figure 23: Parameter fit for tlie intercepts and tlie slopes of tlie curves in Fig. 17 
for the total multiplicities. The numbers over each point in this plot indicate the 
value of the extra-dimensions. The {a, q) parameters are fitted to a straight line 
q = a a + f3 independently of the numbers of extra dimensions, (a) is the fit for 
5,500 m, (b) for 15,000 m. The benchmark is also shown in the plot, but has not 
been used in the fit. 
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Figure 25: Plot of the average radius R of the core of the shower of photons as a 
function of Ei for a black hole with a varying number of extra dimensions. The 
benchmark result is also shown for comparison, (a) 5,500 m, (b) 15,000 m. 
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Figure 26: Plot of the average size of the core R of the shower of e as a function 
of El. (a) 5,500 m, (b) 15,000 m. 
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Figure 27: Plot of the average size R of the core of the shower of // as a function 
of El, (a) 5,500 m, (b) 15,000 m. 
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Figure 28: As above for tt 
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Figure 29: The proton core size as a function of Ei, (a) 5,500 m, (b) 15,000 m. 
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In analogy to Figures 23 and 24, we show in Figure 30 that for a given 
setup there is a hnear relation between slopes and intercepts of Eq. (90) 

q' = a'a' + (3' (91) 

with a' and /9' typical for a given setup, but again independent of n. 



7 Discussion 

A rather detailed analysis was presented of some of the main observables which 
characterize the air showers formed, when a high energy collision in the atmo- 
sphere leads to the formation of a mini black hole. We have decided to focus 
our attention on the particle multiplicities of these events, on the geometrical 
opening of the showers produced, and on the ratio of their electromagnetic to 
hadronic components, as functions of the entire ultra high energy spectrum of 
the incoming primary source. We have shown that in a double logarithmic scale 
the energy vs multiplicity as well as the energy vs shower-size plots are linear, 
characterized by slopes which depend on the number of extra dimensions. We 
have compared these predictions with standard (benchmark) simulations and cor- 
rected for the energy which escaped in the bulk, or emitted by the black holes 
at stages prior to the Schwarzschild phase. Black hole events are characterized 
by faster growing multiplicities for impacts taking place close to the detector; 
impacts at higher altitudes share a similar trend, but less pronounced. The mul- 
tiplicities from the black hole are larger in the lower part of the energy range, 
while they become bigger for higher energies. We should also mention that, given 
the choice made for our benchmark simulations, here we have been considering 
the worst scenario: in a simulation with an impacting neutrino it should be pos- 
sible to discern between the two underlying events, whether they are standard 
or black hole mediated. The lateral distributions appear to be the most striking 
signature of a black hole event. Due to the higher p^s involved, they are much 
larger than in the benchmark standard simulations. 

Our analysis can be easily generalized to more complex geometrical situations, 
where a slanted entry of the original primary is considered and, in particular, to 
horizontal air showers, which are relevant for the detection of neutrino induced 
showers, on which we hope to return in a future work. 

To address issues related to the Centauros, one has to concentrate on black 
holes produced closer to the detectors, since the main interaction in these events is 
believed to have taken place somewhere between and 500 m above the detectors 
[57]. The multiplicities in the showers will of course decrease and are expected, 
based on the above graphs, to get very close to the values observed, while the 
ratio of iVem/A^hadron will also approach the observed values, since the kaons pro- 
duced by the black hole's "democratic" evaporation will not have enough time 
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Figure 30: Parameter fit for tlie curves in Fig. 25, describing tlie openings of the 
showers of photons, (a) is the fit for 5,500 m, (b) for 15,000 m. 
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to decay and will be counted as hadrons. Another crucial question in connection 
with the mini black hole interpretation of the Centauro events is whether the 
partonic decay products of the black hole passed through a high temperature 
quark-gluon plasma phase and whether a Disoriented Chiral Condensate (DCC) 
was formed, before it finally turned into hadrons. In the present analysis it was 
implicitly assumed that no such phase was developed. However, if a DCC forms, 
the prediction for the ratio of the electromagnetic to hadronic component will be 
drastically different. 

Needless to say, the subject of black hole production and evaporation touches 
upon several different subfields of theoretical and experimental high energy physics 
and astrophysics. The fields of Cosmic ray physics and of the air shower forma- 
tion in the atmosphere; of quantum gravity/string theory, of non-perturbative 
low energy QCD and of the quark-gluon plasma phase, just to mention a few, 
and this may not even be a complete list. Given our incomplete understanding of 
all of these, it is clear that a lot more has to be done, before one can safely com- 
pare the theory to the observational data. The analysis of diffractive interactions 
at those energies, in particular, will require considerable attention. Nevertheless, 
in our view, the issues involved are very important and deserve every effort. 
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POSTSCRIPT 

After completing these studies we have been informed by D. Heck that a 
new version of CORSIKA has been released, which is able to deal with neutrino 
primaries. As we have discussed in our work, the use of a more realistic benchmark 
based on neutrino primaries does not invalidate the results of our simulations, 
but should actually render the differences between mini black hole mediated and 
standard events even more pronounced. 
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